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An approach to correlated dynamics of quantum nuclei and electrons both in dynamical interac- 
tion with external environments is presented. This stochastic quantum molecular dynamics rests 
on a theorem that establishes a one-to-one correspondence between the total ensemble-averaged 
current density of interacting nuclei and electrons and a given external vector potential. The theory 
allows for a first-principles description of phenomena previously inaccessible via standard quantum 
molecular dynamics such as electronic and nuclear relaxation in photochemistry, dissipative corre- 
lated electron-ion dynamics in intense laser fields, nuclear dephasing, etc. As a demonstration of the 
approach, we discuss the rotational relaxation of 4-(N,N-dimethylamino)benzonitrile in a uniform 
bath in the limit of classical nuclei. 
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Quantum molecular dynamics (QMD) approaches, as 
formulated, e.g., in the Born-Oppcnhcimer, Ehrenfest or 
Car-Parrinello methods, have proved to be extremely use- 
ful tools to study the dynamics of condensed systems |l|, 
0, 0|- In these approaches, the time-evolution of the nu- 
clear degrees of freedom is described by classical Newton- 
type equations, and the forces acting on the nuclei are 
typically derived from the electronic wave- function at the 
instantaneous nuclear configuration. By construction, 
Born-Oppenheimer and Car-Parinello QMD refer to the 
electronic ground-state energy, whereas Ehrenfest QMD 
also allows to take excited electronic wave-functions into 
account. 

In these approaches, energy dissipation and thermal 
coupling to the environment are usually described with 
additional thermostats coupled directly to the classical 
nuclear degrees of freedom. Alternatively, Langevin dy- 
namics can be employed for the time-evolution of the 
classical nuclei Q. Other routes that take environment 
effects into account are provided by the so-called quan- 
tum mechanics/molecular mechanics (QM/MM) meth- 
ods or by implicit solvation or continuum models (see, 
e -g-j [a])- Here, the system of interest is treated quantum 
mechanically (e.g., with one of the above QMD methods) 
and the environment is treated explicitly within classical 
molecular mechanics or, implicitly, via an electrostatic 
continuum model. 

However, all of the above approaches share the com- 
mon feature that the electronic subsystem is treated al- 
ways as a closed one, namely it cannot exchange en- 
ergy and momentum with the environment (s). This is 
a major limitation since in many physical situations, 
electrons, being lighter particles than nuclei, respond 
much faster than the latter ones to dynamical changes 
in their surrounding. This can lead to electronic transi- 
tions among excited states (manifested, e.g., in an effec- 
tive electron temperature rather different than the ionic 
one Q ) , which in turn may result in different ionic forces 
compared to the forces obtained from a closed electronic 



quantum system at (typically) zero temperature. An 
even more complicated situation may arise (for instance, 
for light nuclei in intense laser fields) when the quantum 
nature of nuclei correlates with the electron dynamics, 
with the external environment (s) mediating transitions 
between electron-ion correlated states. 

In this paper we propose a novel approach that allows 
to treat both the electronic and (in principle, quantum) 
nuclear degrees of freedom, open to one or more environ- 
ments. Our approach, which we term stochastic quantum 
molecular dynamics (SQMD), is based on a stochastic 
current density-functional theory for the combined dy- 
namics of quantum electrons and nuclei. The theorem 
of stochastic time-dependent current-density-functional 
theory for a single particle species has been proved in 
Ref. @, 0. Here we extend it to its non-trivial multi- 
species case. In the limit of classical nuclei, this for- 
mulation reduces to a molecular dynamics scheme which 
couples the nuclear degrees of freedom to an open elec- 
tronic system that can exchange energy and momentum 
with the environment. 

In the following, we consider a system of N e elec- 
trons with coordinates r = {r^} and N n — Yl s Ns,n nu- 
clei, where each species s comprises iV^ particles with 
charges Z s j, masses M SJ , j — l...N Stn , and coordi- 
nates R = {R s .j}, respectively. To simplify the notation 
we denote with x = {R, r} the combined set of elec- 
tronic and nuclear coordinates and we use the combined 
index a = {s, j} for the nuclear species. Now, suppose 
that such a system, subject to a classical electromagnetic 
field, whose vector potential is A(i), is described by the 
many-body Hamiltonian Hs(t) 



H S (t) =T e (r, t) + Wec(r) + £7 e xt,e(r, i) + W ea (r, R) 
+ T n (R, t) + W nn (R) + #ext,n(R, t), 



(1) 



where T (t) and T n (t) are the kinetic energies of elec- 
trons and ions, with velocities Vk(t) = [pk + eA(rk,t)]/m 
and V a (t) = [P a - Z a A(R a ,t)]/M a , respectively. The 
electron-electron, electron-nuclear and nuclear-nuclear 
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interaction terms take the form 
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and we also allow for external time-dependent scalar po- 
tentials £/ext,e(l, t) and ?7cxt,n(IL that act on electronic 
and nuclear degrees of freedom, respectively. We then let 
our system interact with one or more environments repre- 
sented by a dense spectrum of bosonic degrees of freedom, 
described by the Hamiltonian Hb (t) . The system and the 
environments are coupled by a bilinear interaction term 
Hsb = "}ZjSjBj. Here S, B, refer to subsystem and 
bath operators, respectively [l7(. The total Hamiltonian 
H(t) of system and environment then takes the form 



H(t) = H s (t)+H B (t) + '£s j B j . 



(3) 



By projecting out the bath degrees of freedom, one then 
obtains the reduced dynamics of the system of inter- 
est [ij]. While, quite generally, we could work with the re- 
sulting non-Markovian dynamics, for the purpose of this 
paper, we further make a memory-less approximation for 
the bath. This leads us to consider the following stochas- 
tic many-body Schrodinger equation (fi = 1) [l8( 

*&*(x,f) = ffs(t)*(x,t) - -iS^S^(x,t) +l(t)SV{x,t), 

(4) 

where, for simplicity, from now on we consider only 
the coupling S to a single bath, which could be po- 
sition and/or time dependent. The function l(t) de- 
scribes a stochastic process with zero ensemble aver- 
age l(t) = and (^-autocorrelation l(t)l(t') = S(t - t'). 
Here rrr describes the statistical average over all mem- 
bers of an ensemble of identical systems with a com- 
mon initial state ^(x, t = 0) (which need not be pure) 
evolving under Eq. (j4]). Given the set of potentials 
Uext(x,t) = t/ext,e(r,i) + l/ext,n(R,*) and A(t) one can 
always find a gauge transformation A(x, t) so that the 
scalar potentials vanish at all times, and implying that 
A(t) = A(x, t). In the following we assume that such a 
gauge transformation has been performed. 

At this stage we face a large number of choices for the 
definition of densities for the construction of a density- 
functional description of this problem. In fact, certain 



groupings of nuclear particle and current densities might 
be better suited to specific physical situations. To main- 
tain flexibility we do not specialize at this point and we 
prove a theorem for the total current density of electrons 
and nuclei. To that end, we introduce electronic and 
nuclear current operators in terms of the velocity fields 
v k (t) and V a (t) 



j(r,t) 



e 
2m 



(5) 



J a (R,t) = ^- J2 {%(t),6(R--kp)} (6) 

Z a = Z ,M a =Mp 

where {p, q} = pq + qp denotes the usual anticommu- 
tator. Likewise, the usual charge density operators can 
be defined as n(r, t) — J2k — ^k) for the electrons 

and N a (R,t) = Y.p-z^z,, M a =M fi S ( R - for each 
nuclear species. The total particle and current density 
operators of the system can then be written as iV(x, t) = 
*) + E« &a (R, t) and J (x, t) = j (r, t) + £ a J a (R, *) , 
respectively. 

Solutions of the stochastic Schrodinger equation, 
Eq. (|4]), leads to an ensemble of stochastic quantum 
trajectories which have ensemble-averaged total charge 

N(x,t) = (N(x,t)) and current densities J(x, t) — 

( J(x, t) ). Here ( • • • } denotes the quantum mechanical 
average. 

Contrary to the multi-component formulation of 
Kreibich et al. 1^, U| we will not employ a body- fixed 
frame for the solution of Eq. ((4]) . This is motivated by the 
fact that for the initial- value problem of Eq. (j4|) a general 
initial condition for the state vector and, in addition, a 
general time- and space-dependent external vector poten- 
tial A(x, t) would break the translational and rotational 
invariance of the original operator Hs(t) — U cx t(?£,t) in 
Eq. ([1]). For our purposes it will therefore be sufficient 
to work in the laboratory frame. We have now collected 
all ingredients to state the following result. 

Theorem. — For a given bath operator S, many-body 
initial state ^(x, t — 0) and external vector potential 
A(x, t), the dynamics of the stochastic Schrodinger equa- 
tion in Eq. (j4j generates ensemble-averaged total parti- 
cle and current densities N(x,t) and J(x, t). Under rea- 
sonable physical assumptions, any other vector potential 
A'(x,i) (but same initial state and bath operator) that 
leads to the same ensemble-averaged total particle and 
current density, has to coincide, up to a gauge transfor- 
mation, with A(x, t) [l9j . 

Proof. — The strategy for the proof of this statement 
follows the technique introduced in Ref. [13] and also 
employed in Ref. [2]. It is thus sufficient to provide the 
relevant equation of motion and summarize the rationale 
behind the procedure. The relevant equation of motion 
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is the one for the ensemble-averaged current density 

d d d 

m— t j a (r,t) + M a — J a (R, t) = N(x,t)— A(x, t) 



df 



dt 



J(x, t) x (V x A(x, t)) + { ^ oc (x, t)) + { _F cn (x, t) } 



(An(x,*)> +m(g(r,t)} + ^M Q (C? n , Q (R, £)>. 



(7) 



Here we have introduced the electron-electron, electron- 
nuclear and nuclear-nuclear interaction densities 

•?"ee(x, t) = -}^6(r- rj)WjW ee (rj - Vj) + mV • a ee 
•^"en(x, = - / J ^(x - Xj)VjW np (a:i — a?j) + m7 ■ <7 en 
+ ^M a V-<7 ne , (8) 



where .F nn can be obtained from with the replace- 
ments r — > R, w GO — > Wim, mV • <r co — > J2 a M a V 
The stress tensors are 



CC i,J 
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(9) 



with a nci j and (T nnij obtained by the exchange v <-> V. 
The remaining terms in Eq. ((7]) describe the electromag- 
netic contribution (terms containing A) and the force 
densities due to the bath 



Q(r,t) = S^j(v,t)S -^S,j(r,t)} 
g a {R,t) = &J a (R,t)S- ^{&S, J(R,i)}. 



(10) 



We can now consider another (primed) system with 
different initial condition, interaction potentials w' ccl w' nn , 
w' en and external vector potential A'(x, i). By taking the 
difference of the equations of motion for the total current 
densities J(x, t) and J'(x, t) in the unprimed and primed 
system we arrive at an equation of motion for AA(x, t) = 
A(x, t) — A'(x, t). The difference of the vector potentials 
AA(x, t) can then be shown to be completely determined 
by the initial condition and a series expansion in time 
about t = 0. If the two total current densities of the 
primed and unprimed system coincide then the unique 
solution - up to a gauge transformation - is given by 
AA(x, t) = 0, when the initial conditions and interaction 
potentials are the same in the two systems, which proves 
the theorem. 

Discussion. — With this theorem we can now set up a 
Kohn-Sham (KS) scheme of SQMD where an exchange- 
correlation (xc) vector potential A xc (functional of the 



initial state, bath operator, and ensemble-averaged cur- 
rent density) acting on non-interacting species, allows to 
reproduce the exact ensemble-averaged density and cur- 
rent densities of the original interacting many-body sys- 
tem. The resulting charge and current densities would 
thus contain all possible correlations in the system - if we 
knew the exact functional. As mentioned before, alterna- 
tive schemes could be constructed (based on correspond- 
ing theorems that could be proved as the above one) 
by defining different densities and current densities (e.g ., 
by lumping all nuclear densities into one quantity [10(). 
While this may seem a drawback, it is in fact an advan- 
tage: some of those schemes may be more appropriate for 
specific physical problems. Irrespective, one would need 
to construct xc functionals for the chosen scheme. While 
this program is possible, it is beyond the scope of the 
present paper. Therefore, as a first practical implemen- 
tation of the proposed SQMD, we consider the limit of 
classical nuclei. This is definitely the simplest case, but 
it is by no means trivial, and by itself already shows the 
wealth of physical problems one can tackle with SQMD 
and the physics one can can extract from it. 

In the limit of classical nuclei we then solve the stochas- 
tic time-dependent KS equations with given KS Hamil- 
tonian HKs{?£,t) for the electronic degrees of freedom, 
at each instantaneous set of nuclear coordinates. The 
result provides the KS Slater determinant ^ifs(x, t) 
from which the forces on nuclei are computed as F Q = 
-(^Ks(x,t)\V a H K s{x,t)\$>Ks(x,t)), for each realiza- 
tion of the stochastic process. One then has direct ac- 
cess to both the average dynamics as well as its distribu- 
tion [2(| . As application we consider the rotational relax- 
ation of 4-(N,N-dimethylamino)benzonitrile (DMABN) 
in a uniform bath. The relaxation rate that enters the 
definition of the operator S can be derived in principle 
from the system-bath interaction term Hsb in Eq. Q [9] . 
In the case of DMABN, relaxation rates are available 
from experiment [lq , so that we choose r = 100 fs for our 
simulation which is within the experimentally observed 
magnitude. We choose as initial state for the SQMD 
simulation a rotated dimethyl side-group (S = 15°) of 
DMABN. In Fig. [T] we show the average angle and the 
binned angle distribution of the dihedral angle S for 50 
stochastic realizations for bath temperatures of OK and 
300K as a function of time and compare with the closed 
system solution. The angle distributions in the lower 
panel of Fig. Q] clearly show the temperature-dependent 
relaxation of the rotational motion in the open system 
case. Most importantly, from the transient dynamics of 
the angle distributions, it is also apparent that the system 
is not relaxing uniformly to the equilibrium configuration 
5 = 0. Rather, it approaches equilibrium via a series of 
quasi-bimodal distributions, with the higher temperature 
"smoothing" these distributions. We also emphasize that 
the damping of the nuclear motion originates in our simu- 
lation exclusively from the forces that are calculated from 
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FIG. 1: (Color online) Upper panel: schematic illustration of 
4-(N,N-dimethylamino)benzonitrile and the dihedral angle S. 
Center panel: comparison between the dihedral angle S as a 
function of time for a closed quantum system (dashed line) at 
K and an open quantum system (solid line) at 300 K. Lower 
panel: Distributions of the dihedral angle 5 for 50 members of 
the stochastic ensemble for bath temperatures of K and 300 
K at different times (as indicated by the dotted connection 
lines between the center and the lower panel). 



the electronic open quantum system wave- functions. No 
additional friction term has been added to the nuclear 
equation of motion. It would be thus interesting to ver- 
ify such a prediction with available experimental capa- 
bilities. 

In summary, we have presented a novel quantum 
molecular dynamics approach which we term stochas- 
tic quantum molecular dynamics (SQMD), based on a 
multi-species theorem of DFT for open quantum systems. 
SQMD allows to treat both the electronic and nuclear 
degrees of freedom open to environments, and, in princi- 
ple, it provides all possible dynamical correlations in the 
system. In particular, SQMD takes into account energy 
relaxation and dephasing of the electronic subsystem, a 
feature lacking in any "standard" MD approach. This 
opens up the possibility to study a wealth of new phe- 
nomena such as local ionic and electronic heating in laser 
fields, relaxation processes in photochemistry, etc. 
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Note that both the electron-phonon and the electron- 
photon interaction can be cast in the form of Hsb- 
The reasons why we work with a stochastic Schrodinger 
equation and not with a density-matrix approach are due 
to both the possible loss of positivity of the density ma- 
trix during dynamics, and the fact that a Kohn-Sham 
Hamiltonian depends on the density and/or current den- 
sity, and is thus generally stochastic 0, |g|. 
As in [7|, |8| we are implicitly assuming that given an 
initial condition, bath operator, and ensemble-averaged 
current density, a unique ensemble-averaged density can 
be obtained from its equation of motion. Therefore, the 
density is not independent of the current density, and our 
theorem establishes a one-to-one correspondence between 
current density and vector potential. 
We have implemented a quantum-jump algorithm [Jj| 
in the real-space real-time package octopus [14| to in- 
tegrate the KS equations of motion. We have used the 
bath operators as in Ref. [15j, and the adiabatic local- 
density approximation for the xc potential. We note that 
the jumps introduced by the bath in the quantum-jump 
algorithm appear like a "surface hopping". The relax- 
ation rates derived from Hsb and the jumps induced by 
the bath, however, provide a solid framework for a sur- 
face hopping scheme. More details of the numerics will 
be given in a forthcoming publication. 



